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Transportation networks are inevitably selected with reference to their global cost which depends 
on the strengths and the distribution of the embedded currents. We prove that optimal current 
distributions for a uniformly injected d-dimensional network exhibit robust scale-invariance prop- 
erties, independently of the particular cost function considered, as long as it is convex. We find 
that, in the limit of large currents, the distribution decays as a power law with an exponent equal 
to {2d — — 1). The current distribution can be exactly calculated in d = 2 for all values of 

the current. Numerical simulations further suggest that the scaling properties remain unchanged 
for both random injections and by randomizing the convex cost functions. 

PACS numbers: 89.75.Hc, 89.75.Da, 89.75.Kd, 89.75.Fb, 05.65.+b, 45.70.Vn, 68.70.+W 



I. INTRODUCTION 

Finding efficient ways of distributing (or collecting) 
matter injected through a given region, spanned, e.g., 
by a regular lattice, from (at) a unique source (sink) is 
relevant to a variety of problems arising both for natural 
and artificial systems. The main mechanisms that arc 
known to achieve a capillary distribution arc by diffu- 
sion or through a network-like structure providing a near 
uniform spatial supply (drainage), or a combination of 
them. Network arrangements, of this kind, are observed 
in many living organisms, like for instance circulatory 
and lymphatic systems in animals or xylem and roots in 
vascular plants, and are also widely employed in artificial 
systems such as electrical or hydraulic transmissions and 
fiuvial basins 0,0,11,0. 

Certainly, one wonders what is the basic selection prin- 
ciple that favors, say, tree- like versus looping network 
structures, in view of the widespread occurrence of both 
forms in nature and elsewhere [1, Q 0- To that end, 
it has been previously shown that tree-like structures 
emerge as local minima of global energy expenditure in 
networks whose transportation cost is physically con- 
strained to be a concave function, such as in the case 
of river networks Q . 

Transportation costs generally depend on the strength 
of currents and on network topology. The best known ex- 
ample is the electrical resistor network 0, [l^ : consider a 
square lattice where a resistor is placed at every bond be- 
tween each pair of nearest-neighbor nodes. Each node is 
externally supplied by a unit flux. If a sink collects all the 
currents, one is capable of controlling the current fluxes in 
all bonds by assigning the potential differences between 
pairs of nodes. The total cost in the transportation - the 



dissipated power - is proportional by Ohm's law to the 
sum of the square of currents: £ {{lb}) = J2b l^fcl^- ^P" 
timal resistor networks are thus obtained by setting the 
potentials in order to minimize the total transportation 
cost. The most immediate way to generalize the resistor 
network case is to replace the exponent "2" by a generic 
exponent 7 > 0. 

The case 7 < 1 has been characterized exactly 0, H, [HI, 
m, • It was shown there that the cost function admits 
many local minima corresponding to configurations with 
currents present only on the bonds of spanning trees 
The case 7 = 1/2 exhibits scaling behavior akin to river 
networks 0, [HI, [H, [II ■ 

The case 7 = 1 is related to the Voter model [l^, mass 
aggregation [H, [13, [Hj [3 , directed sandpilc models [20t . 
and Kleiber's law of metabolic scaling of living organism 
[2ll |. Recently optimal transportation networks with a 
global constraint have been studied in a variety of con- 
texts [22I, [2^ , where different topologies of the optimal 
transportation network arising in the cases 7 < 1 and 
7 > 1 have been investigated also by Bohm and Mag- 
nasco [23} . 

The case 7 > 1 is an example of convex transportation 
costs. Important real-world examples are found, for ex- 
ample, in road traffic analyses where more cars cause dis- 
proportionately higher costs (travel times) , usually mod- 
eled as convex functions or in any electricity distribu- 
tion network owing to Ohm's law. Convex cost functions 
also occur in many operations research applications as 
pointed out in [2^ . 



II. NUMERICAL RESULTS 

Here we address the current distribution for the case 
corresponding to a general convex cost function of the 
circulating currents: 
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found to be: 



FIG. 1: The current conservation at node x. The sum of 
aU currents flowing into node x must be equal to the sum of 
currents flowing out of node x: Ib^ + Ib2 — hs — hi + ix = 0- 



where the sum spans all network bonds b. In specific 
examples and numerical simulations we will consider the 
particular class of convex functionals with E{\ |) =| 
lb /7 and 7 > 1. It will be shown, however, that 
our findings are valid for an arbitrary convex cost func- 
tion E with finite first derivative which depends only on 
current strength. Notice that for a non-linear resistor 
network the transportation cost £ is proportional to the 
dissipated power only for the case E{\ If, |) oc| /& as 
discussed in [2^,[2^ . We will focus here only on finite net- 
works. The goal is to determine the current probability 
distribution corresponding to the current configuration 
{/;,} which minimizes the total transportation cost £ in 
the large size limit and its scaling behavior as a function 
of the cost E{\ I I). 

The minimization of £ is subject to the local constraints 
of current conservation at each node x: the sum of all cur- 
rents flowing into a node, taken as negative (positive) if 
directed outward (inward), must equal the injected nodal 
flux, ix (see Fig. [T]). 

Because such constraint is linear in the /f,s, if iS is a 
convex function then it admits a single global minimum. 
The existence and uniqueness of the solution for infinite 
networks have been addressed elsewhere [2§ |. 
As shown below, we find quite generally that the cumu- 
lative probability distribution function (CPDF), i.e., the 
fraction of currents |/{,| larger than /, obeys the finite- 
size scaling 



Vd,7 > 1, / > 



(2) 



where L is the linear size of the system and F is the 
scaling function which at large x behaves as F{x) ~ x^~'^ 
with T = {2d - l)/{d - 1) and lim^^o ^"(2;) = 1. This 
scaling form holds independently of the particular convex 
cost function considered. The standard scaling behavior 
for the case E{\ /f, |) =| Ih \'' h and 7 < 1 [ll|, 113 
corresponding to a non-convex cost function, was instead 
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7 < 1, ^ ^ 1 



(3) 



with lim^-to /(a;) = const, i.e., a pure power law is ob- 
tained in the large size limit ( r = 1.43 ± 0.03 in d = 2 
and 7 — 1/2). The 7=1 case was solved exactly in all 
dimensions using a mapping to reaction diffusion models 
[l6t and r = 2(d -I- l)/((i -I- 2) when the dimensionality 
is lower than the upper critical dimension, dc = 2, and 
T = 3/2 when d > 2. In the present case no upper criti- 
cal dimension is found above which the exponent remains 
the same. 

Let us first describe the results of the numerical determi- 
nation of current configuration minimizing Eq.lll) with 
E{\ lb I) =1 lb I'' and 7 > 1 in a square lattice (i.e., 
d = 2) of linear size L where a uniform input at each 
site = 1 is assumed (at the sink where all currents 
are collected one has isink = —L'^+1). We have used 
open boundary conditions for simplicity because we do 
not expect that they influence the scaling behavior in 
the large size limit. Because our problem reduces to the 
minimization of a convex function of many variables, we 
have used the nonlinear conjugate gradient method [28| . 
The CPDF, PciI\L), is plotted in Fig. [H for lattices of 
different sizes L and for 7 = 2, the resistor network. It 
has the following scaling behavior: 
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for < / < L 
for L < / < 



(4) 



with T = 2.975 ± 0.045. In the inset of Fig. [Uthe CPDF 
is plotted versus 1/ L for various L, showing that indeed 
Pc{I\L) is a homogeneous function of the ratio I/L. 

We have also performed numerical optimizations with 
different values of 7 > 1 in c? = 2. Figure [3] shows 
pictures of the optimized current configuration and the 
corresponding CPDFs are plotted in Fig. [D Because 
data overlap, and although the current configuration of 
global minimum for £ varies with 7 as shown in Fig. [3l 
it is suggested that the distribution of currents is inde- 
pendent of the exponent 7 when 7 > 1. Additionally, we 
performed simulations on a 2 — d triangular lattice. The 
scaling behavior of the CPDF proved to be independent 
of the underlying lattice's structure. 



III. ANALYTICAL RESULTS 

We attack the problem analytically in the continuum 
limit (this will be justified a posteriori). We begin to il- 
lustrate the procedure in detail for the case E{\ lb |) =| 
lb I'' /7; later on we will extend it to the more general 
case. 

Under the assumption that the current distribution, in 
the large size limit, does not depend on the shape of 
the volume we enclose our system in a region n = 
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FIG. 2: (Color online) CPDF for square lattices of different 
sizes (L — 201 (red circles), 301 (green squares), 501 (orange 
rhombi), 701 (blue triangles)) and for 7 = 2 (on log-log scale). 
Fitting the tail of the distribution with a power law yields a 
value for the exponent r ~ 3, e.g. for L = 301 we get t — 
2.975 ± 0.045 (the exponent of the power law was estimated 
using the method of maximum likelihood [29l . |30||. and the 
error was calculated with the bootstrap method ^l[). In the 
inset we plot P{I\L) vs I/L, for all L's and J's considered 
above. The collapse of the curves indicates that the CPDF is 
a function of the ratio I/L. 



{x : |lx|l^ < L}, whose volume will be denoted as 
We have defined the norm as ||x||^ = jx^l'^]-^/''. 
In the case 7 = 2, 17 is a sphere of radius R = L. 
Let j(x) be the current density at location x whereas 
i(x) = — |r2|i5'' (x)] is the external input. The Dirac 
delta distribution 6'^ (x) [a compact notation for the d- 
dimensional notation 5"^ (x) = 6{xi) . . . 6{xd) ] represents 
the sink at the origin, whereas ia is the uniform input. 
The components j^(x) {fi = 1, d) of the vector j(x) 
represent the currents along the positive direction of co- 
ordinate axes at position x. 

In the continuum case we define a cost functional anal- 
ogous to Eq. Il]) as: 



d X 



7 



(5) 



We search for the current configuration that minimizes 
the cost function Eq.([5]) with the constraint of the cur- 
rent conservation law, V ■ j(x) = i(x), at each position, 
X. This is done by introducing a Lagrange multiplier 
(potential), l^(x), at each position, x, and solving the 
following equation: 



= 



'5j,<(x) 



d'^.T y V • j = 



b-^(x)|2-7 dx. 



(6) 

where /i = 1, 2, . . . , d. Because we expect that the CPDF 
does not depend on boundary conditions in the large size 



FIG. 3: (Color online) Currents' intensity in the optimal con- 
figuration for different values of the exponent 7 and for the 
size L = 151 of a square lattice. From top left to bottom 
right: 7 = 1.5, 7 = 2, 7 = 4 and 7 = 6. The petal-like 
arrangement of currents is not an artifact of the underlying 
lattice geometry but arises from the decomposition of current 
vectors into components. The direction of currents is towards 
the center (directed networks). The colors indicate the in- 
tensity of currents: Yellow: L < I, Purple: L/2 < I < L, 
Green: L/4 < / < L/2, Blue: L/8 < I < L/4, Black: 
L/16 < / < L/8, Red: / < L/16. 



limit, we choose j(x) = at the boundary. We now as- 
sume that the solution depends only on ||x||^ : because 
£ is convex, a solution with this property (if it exists) is 
the unique solution. Using the above choice of the in- 
put currents d'^x i(x) = and the conservation law 
V • j(x) — i(x), by applying Gauss theorem one gets 
that the boundary condition on the volume are auto- 
matically satisfied. Equation ^ together with current 
conservation gives the following radial current density as 
the optimal solution of Eq. ([5]) 



J(x)=x- 



1 - 




(7) 



This solution is radially symmetric only with the met- 
ric defined in terms of the 7-norm itself: thus the equi- 
currents lines defined by ||j||-y = constant are circles 
only for 7 = 2, whereas when 7 —>■ 00 (1) they be- 
come squares with sides parallel to the coordinate axis 
(45°-tilted squares). The CPDF is given by Pc{j\L) = 



d'^x 



(bM('^)l ~ j) if consider currents' 
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FIG. 4: (Color online) Comparison of CPDFs of global min- 
imum configurations for the class of convex functionals ([T]) 
with E{\ Ifj I) =1 P /7 and different 7-values: 7 = 1.5 (red 
circles), 2 (blue squares), 4 (green rhombi), 6 (orange trian- 
gles). The black solid line is eq.(|9}, the CPDF of the analytic 
solution of eq.®. The shapes of the network boundaries in 
the various cases are chosen as explained in the text. 



components, or by 

PML)^m-' I d'^xe(|jj(x)||,-j) (8) 

Jo 

if we consider current norm [s^l [0(2;) = 1 if 2; > and 
zero otherwise]. It can be shown that asymptotic be- 
haviors arc the same in both cases. Using the explicit 
solution Eq.([71) one sees that Pc{i\L) depends only on 
the dimensionless ratio j / (ioL) for all d. When d = 2 
Eq.([5]) takes the simple closed form: 

Pcim^pl^^y Fiz) = [VT+^ - z)' (9) 

which is of the kind anticipated in Eqs.ijl]) and ([2]) with 
r = 3. The prediction Eq.([ni) is shown in Fig. d] com- 
pared to CPDFs of numerical simulations calculated con- 
sidering currents' components. Even though we do not 
have the explicit analytical form for d > 2 it is not dif- 
ficult to verify the asymptotic behavior of Eq.lU). In- 
deed for j/L 1 the leading contribution in Eq.® 
comes from ||x/L||^ <C 1 leading to Pc{j\L) (j/i)*^"^) 
with T = {2d — l)/{d — 1) for d > 1. The minimum 
||x||^ is given by the underlying lattice spacing, say a, 
and so, according to Eq.([7]), the maximum current is 
of order aio{L/a)'^. Thus scaling holds in the region 
1 < j/ioL < [L/aY^^. The special case d = 1 is triv- 
ial and one gets Pc{j\L) = 1 - 2j/{iQL)Q[l - 2j/{ioL)). 
Notice that in this case the scaling region has shrunk to 
zero. 

The case 7 < 1 cannot be treated in the same way as 
above. In fact in it was shown that for the functional 



([T]); with E{z) oc z'', any current configuration, with cur- 
rents being different from zero only on the bonds of a 
spanning tree, is a local minimum. For such solutions 
the second term in Eq.® would diverge in correspon- 
dence of the bonds not belonging to the spanning tree. 
The above results can be generalized to the case of a 
generic convex cost function with finite first derivative as 
follows. 

While in the presence of an underlying network the 
natural choice for the cost function is given by eq. ([1]) in 
the continuum there are at least two natural choices. The 
first choice corresponds to £ = J^-^d'^x X^^^dj^l^)!) 
whereas the second one is given by: 

£ = f d'x Ei\\j{^)\\,). (10) 

It can be shown that, although these two functionals have 
different optimal configurations, their CPDF have the 
same scaling behavior for small and large currents. For 
the previous case, E{z) = z'' /'j, the two choices coin- 
cide. The minimum of the cost function, Eq. ljlOp . in the 
domain fl = {x : ||x||-), < L} and with the constraint of 
current conservation proceeds as before. The stationarity 
conditions of the constrained problem are: 

E' (||j(x)||,) M^j,(x) = V (||x||,) If^x, 

(11) 

where, as before, we have assumed that the potential 
V{x.) is a function of l|xl|^. The solution is given by 
j^(x) = a;p/(||x||^) with V{z) satisfying the equation 
E'{f) = V'{z) (prime indicates the derivative with 
respect to the argument). Imposing current conservation 
we get f{z) = io/d{l — [L/ zY) leading again to solution 
(l7|). In turn this implies that the scaling behavior of the 
CPDF as defined in Eq. ([S]) is independent of the specific 
cost function as long as it remains convex. 



IV. INHOMOGENEOUS CASES 

We have further tested the robustness of our results by 
performing additional numerical simulations on systems 
subject to independent, equally distributed random cur- 
rent injection, z(x) > at the nodes, or in the presence 
of non-uniform conductivity where the cost function is 
given by £({/(,}) = X^h^fcl^fcP where fcf, are random 
positive numbers. 

As random distribution for injections and conductances 
we have chosen a power law to ensure a high degree of in- 
homogeneity. The simulation results of Fig. [5] show that 
the leading trend for large currents remains the same as 
in the uniform case studied above. Thus it is plausible 
that the scaling behavior of the current distribution cor- 
responding to the optimal solution of the uniform case 
might remain the same even for the more general case of 
a spatially varying convex cost function. These results 
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FIG. 5: (Color online) Comparison among the CPDF for the 
uniform case (red circles) with two examples of the heteroge- 
neous conductivity (green rhombi) and injection cases (blue 
triangles) for a L — 151 lattice, with E{\ h \) =\ h \~' /7 and 
7 = 2. The probability distribution used to extract random 
resistances and injections is a power law with exponent —1.5, 
at large values, in order to provide a high degree of inho- 
mogeneity. Ic is properly chosen to show that the power law 
exponent is the same at large currents for the three configura- 
tions. The shown straight line has slope —2 as our analytical 
results predict. 



differ from the case of random transportation dynamics 
[13, [3 fo'^ which it was shown that the uniform injec- 
tion case is equivalent to our optimization problem with 



7 = 1 [HI- Indeed for these models the scaling behavior 
of the CPDF proves sensitive to the distribution of the 
injections. However the present numerical results sug- 
gest that this equivalence can not be generahzed to the 
random injection case. 



V. CONCLUSIONS 

In summary, we have studied a class of optimal 
transportation networks with a convex cost function as 
given by Eq.([T]) whose prototype is £ {{h}) = X^b I P 
with 7 > 1. The optimal current configurations exhibit 
a probability distribution function characterized by a 
scaling behavior given by Eqs.Q and The scaling 
exponent of the current distribution proves robust with 
respect to: (i) the choice of the transportation cost, 
as far as it is convex and has finite first derivatives 
with respect to the currents; {ii) the distribution of 
injected currents; [Hi) position-dependent (convex) cost 
functions. The analytical results show that the exponent 
of the asymptotic power-law behavior of the current 
probability distribution function varies continuously 
from 3 in two dimensions to 2 at infinite dimensions 
with no evidence of an upper critical dimension. 
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